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Abstract 

We describe an algorithm that computes the ground state energy and correlation functions 
for 2-local Hamiltonians in which interactions between qubits are weak compared to single-qubit 
terms. The running time of the algorithm is polynomial in n and 6"^, where n is the number 
of qubits, and 6 is the required precision. Specifically, we consider Hamiltonians of the form 
H = Hq + eV, where Hq describes non-interacting qubits, F is a perturbation that involves 
arbitrary two-qubit interactions on a graph of bounded degree, and e is a small parameter. 
The algorithm works if |e| is below a certain threshold value eo that depends only upon the 
spectral gap of Hq, the maximal degree of the graph, and the maximal norm of the two-qubit 
interactions. The main technical ingredient of the algorithm is a generalized Kirkwood-Thomas 
ansatz for the ground state. The parameters of the ansatz are computed using perturbative 
expansions in powers of e. Our algorithm is closely related to the coupled cluster method used 
in quantum chemistry. 
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1 Introduction and summary of results 



Perturbation theory provides a systematic way of getting approximations to eigenvalues and eigen- 
vectors for a variety of quantum spin models. Arguably, a significant part of analytical and numerical 
results of condensed matter physics has been obtained using perturbative expansions in some small 
parameter. Quite recently the methods of the perturbation theory have been successfully employed 
in quantum complexity theory. In Ref. [1] Kempe, Kitaev, and Regev used perturbative reductions 
to show that the problem of computing the ground state energy of a Hamiltonian with two-qubit 
interactions is QMA-complete. After that Terhal and Oliveira |2] generalized this result to local 
Hamiltonians on a 2D lattice. 

Our main goal is to examine whether the methods of the perturbation theory provide an efficient 
computational algorithm for the simulation of quantum spin systems. In this paper we focus on 
the simulation of low-temperature properties, namely computing the ground state energy and spin- 
spin correlation functions for the ground state. An efficient algorithm must have a running time 
T = 0(n° 6~^), where n is the number of spins, 5 is a precision up to which we need to compute the 
ground state energy or a correlation function, and a, (3 > are some constants. 

Before stating the results, let us describe the spin models that we shall consider. Let Q = {C, £) 
be a graph with a set of vertices £, |£| = n, and set of edges £. Suppose n spins-1/2 (qubits) are 
located at vertices u & C and spin-spin interactions are located on edges (m, v) G S. The Hamiltonian 
is 

H{e) = Ho + eV, ifo = A„ V = ^->- (1) 

u€C {u,v)g£ 

Here Vu^v is an arbitrary operator acting on a pair of qubits u, v, and e is a real number. The operators 
Hq and V are called the unperturbed Hamiltonian and the perturbation. We shall always assume that 
Au > for all M G C Accordingly, the unperturbed Hamiltonian Hq has a non-degenerate ground 
state 

= |o,o,...,o), Ho\n) = o. 

Most of the time, all we will need to know about Hq and V are the following parameters 

A = minA„, J= max ||V"u,t,||. (2) 

ueC {u,v)&£ 

The parameter A is the gap between the smallest and the second smallest eigenvalue of Hq, while the 
parameter J characterizes a strength of the perturbation V. Let d be the maximum vertex degree 
of the graph Q, 

d = max I {t> : {u, v) E £} \ . (3) 

The quantity we are interested in is the smallest eigenvalue of H{e), which we shall denote by E{e). 
Clearly, E{e) is a continuous concave function of e and -E'(O) = 0. Besides, since we assume that 
A > 0, the standard perturbation theory arguments [3] show that E{e) is analytic at e = and the 
Taylor series 

oo 

E{e)=Y,E,eP (4) 
p=i 

converges absolutely for ||eV^|| < A/2. The following theorem proved by Yarotsky [1] asserts that 
E{e) is non-degenerate for sufficiently small e and sets a lower bound on the spectral gap. 
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Theorem 1. Suppose |e| < 2eo, where 



Then the smallest eigenvalue E{e) has multiplicity 1 and the gap between E{t) and the second smallest 
eigenvalue of H{e) is at least A/2. 

(The exphcit value of eo has not been stated in Ref. |4j.) We shaU provide an alternative proof of 
Theorem [1] in Sections [2|3l 

As was shown by Osborne in Ref. [5], Theorem [T] implies that expectation values of local observ- 
ables on the ground state of H{e) can be efficiently computed within any constant precision S by 
simulating quantum adiabatic evolution along the path connecting H{0) and H{e). However, the 
running time of such simulation scales exponentially as a function of 6^^. As was noted in Ref. 0, 
it means that simulation of the adiabatic evolution does not yield a polynomial-time algorithm for 
computing the ground state energy. 

The perturbation theory provides an approximation to the ground state energy by truncating the 
series Eq. at sufficiently high order p. In order to understand whether this approach can be used 
to construct an efficient computational algorithm, two separate issues have to be addressed: 

Ql: What is the convergence radius of the perturbative series? 

Q2: What is the computational cost of finding the coefficients in the perturbative series? 

Note that the radius of convergence of the series Eq. (jl]) is a property of the Hamiltonian H{e) only. 
It does not depend upon what particular perturbative expansion has been used to find the coefficients 
Ep. The following theorem allows one to answer the first question. 

Theorem 2. The Taylor series E{e) = Yl'^=i^p^^ converges absolutely for |e| < 2eo. Furthermore, 



Thus if one needs to compute E{e) with a specified precision 6, it suffices to compute the coeffi- 
cients El, . . . , Ep, where p = log2 {n6~^) -|- 0(1) (assuming that A is a constant that does not depend 
on n). 

Answering the second question has nothing to do with the convergence radius of the series Eq. (jlj) 
(as long as it is non-zero). One can compute the coefficients Ep by choosing e so small that ||e V^|| <ti A. 
In this regime the standard perturbation theory is applicable, for example, the self-energy opera- 
tor formalism, see Refs. [H E], or the Rayleigh-Schrodinger expansion, see Ref. [7|. Clearly, the 
computational cost of finding the coefficients Ep varies for different methods. 

In the present paper we compute the coefficients Ep using the Kirkwood- Thomas ansatz for the 
ground state. It was originally proposed in Ref. |8] for translation-invariant Ising-like Hamiltonians 
with a transverse magnetic field. The translation-invariance constraint has been removed in the 
later work by Datta and Kennedy [9j. We use the generalized Kirkwood-Thomas ansatz proposed 
by Yarotsky [1| which is applicable to any spin Hamiltonian with sufficiently weak interactions. It 
allows us to prove the following. 




4 



Theorem 3. Suppose d is a fixed constant independent of n. Then there exists an algorithm with a 
running time riexp {p{p)) that takes as input a triple {Hq, V,p) and outputs Ei, . . . , Ep. 

An immediate consequence of Theorems [2][3] is 

Corollary 1. Suppose A, J, d are fixed constants independent of n and \e\ < eo- Then there exists 
an algorithm with a running time poly{n, 6) that computes E{e) with an absolute error at most 6. 

Besides, it follows from Theorems [2|3] that the energy density E{e)/n can be computed with a 
precision 5 in a time n ■ poly{6~^). 

Note that while computing the coefficients Ei, . . . , Ep we cannot afford the running time to grow 
faster than exp {0{p)) (for fixed n) since we need p ~ log {n6~^) to achieve the desired accuracy. The 
perturbative expansion based on the Kirkwood-Thomas ansatz has two special features that make 
the scaling exp (0(p)) possible: (i) The parameters of the ansatz are complex amplitudes C(M) 
assigned to subsets of vertices M C £. The recursive equations specifying the amplitudes C{M) 
are described by a polynomial of a constant degree, see Section I3.lt (ii) The perturbative expansion 
C(M) = C*p(M) has a property known as the linked cluster theorem, namely, Cp{M) = 

unless M can be covered by a connected subgraph of size 0{p), see Section I^?T1 The number of such 
subgraphs grows only exponentially with p, see Section 14. 2[ It implies that the number of non-zero 
coefficients Cp{M) grows as nexp (0(p)), see Section ISTTl 

Naturally, one could run the algorithm from Theorem [3] to compute the truncated series for E{e) 
even if |e| > eo. The running time will be polynomial in n and as long as |e| is smaller than the 
convergence radius R of the series Eq. (jlj). Although we believe that R must be close to A/{dJ) 
(see a discussion at Section [6]), its exact value cannot be easily found. In practical simulations, one 
could evaluate R by computing sufficiently many coefficients Ep and using the fact that R~^ is the 
largest accumulation point of a sequence \Ep\^/P, p = 1, . . . , oo, see Ref. jTUj. Note that in general 
the singular point (points) of E{e) with |e| = R does not lie on the real axis and thus cannot be 
identified with a quantum phase transition point of H{e) (since we consider finite systems, the latter 
is not even well defined). 

Obviously, efficient computation of E{e) is possible due to the presence of a small parameter e in 
the problem. However it should be emphasized that the condition |e| < eo does not imply that the 
ground state \ip) of H{e) is close to the ground state \Q) of the unperturbed Hamiltonian Hq. In fact, 
one should expect that \ip) and are almost orthogonal for large r^. To illustrate this statement, 
consider as an example the perturbation V = —Jj^uec-^^^ where X is the Pauli operator, and 
the unperturbed Hamiltonian Hq = A ^^^^ Clearly, the ground state of H{e) is a product 

of one-qubit states, = <S>u€c\'^u)- ^ simple calculation shows that {0\ipu) = cos (6*72), where 
cos (6) = (l + 4e^ J^/A^)~^/^. Thus for any fixed e the overlap = (cos {0/2)Y gets exponentially 

small as n increases. However the reduced density matrices of the ground states and for any 
subset of qubits of constant size are indeed close to each other for small e. In other words, for small 
e the state IV') describes small density quantum fluctuations of the background state \VL). One could 
speculate that this statement remains true for arbitrary weak perturbations as well. The Kirkwood- 
Thomas ansatz for the ground state of H{e) used in the present paper provides a convenient way 
to quantify the "density of quantum fluctuations" and prove that it is indeed small for |e| < eo. 
Unfortunately, our approach does not allow us to make any statements about the validity of the area 

-'^This effect is analogous to the well-known "orthogonality catastrophe" observed by Anderson in Ref. [TT] for 
non-interacting fermions in a presence of a scattering potential. 
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law or to decide whether the ground state can be well approximated using the PEPS ansatz, see 
Ref. [12] . 

The Kirkwood-Thomas ansatz is not well suited for computing spin-spin correlation functions 
because it provides an unnormalized ground state. We avoid this problem using the standard relation 
between the correlation functions and the linear response of the ground state energy to a small 
perturbation. It allows us to prove 

Theorem 4. Let Ou,v be a Hermitian operator acting non-trivially only on qubits u,v & C. Suppose 
\\Ou,v\\ ^ 1- The expectation value of Ou,v on the ground state of H{e) can be computed with a 
precision 6 in a time T = poly{6'^) as long as \e\ < eo/2{d+ 1). 

Remark: After the present work has been completed, it was communicated to us by F. Verstraete [I3] 
that the simulation algorithm based on the Kirkwood-Thomas ansatz is closely related to the coupled 
cluster method originally introduced by Coester [H]. The coupled cluster method is extensively used 
for numerical simulations in quantum chemistry, see a review [15], as well as in condensed matter 
physics, see a review [16] and the references therein. Accordingly, from the perspective of practical 
simulations, the algorithm described in the present paper is certainly not a new one. However, we 
believe that our results provide the first rigorous proof that the coupled cluster method yields a 
polynomial-time simulation algorithm for spin Hamiltonians with weak interactions. 

The rest of the paper is organized as follows. Section [2] provides the necessary background on 
the generalized Kirkwood-Thomas ansatz. It mostly follows Ref. [4], although some of our proofs are 
technically different (in particular. Lemma Hj). Section [3] shows how to solve the Kirkwood-Thomas 
equations using a power series and proves Theorem O In Section H] we prove that our perturbative 
expansion obeys the well-known linked cluster theorem and establish an upper bound on the number 
of linked clusters on a graph. The algorithms for computing the ground state energy and spin-spin 
correlation functions are explicitly described in Section |5] which provides a proof of Theorems [3|^ 
Some open problems are discussed in Section O Appendix A contains a technical lemma proving 
submultiplicativity of the norm of creation operators. 

2 Kirkwood-Thomas ansatz for the ground state 
2.1 Creation operators 

Define one-qubit operator = |1) (0|. Let a]^ be the operator on qubit u tensored with the identity 
on all other qubits. For any non-empty subset of vertices M <Z C denote aj^^ = j^^gj^^ aJi. Note that 
the operators aj^^ are nilpotent, (aj^^)^ = 0, and that they pairwise commute: ajy^aj^ = a^^ajy^. Also, 
one can easily check that the operators {ajy^}, 7^ M C £ are linearly independent. (All these 
definitions and properties apply to a and gm operators as well). 

Definition 1. A creation operator is an operator that can be written as 

C= J2 CiM)al 

for some complex numbers C{M). 
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For any given creation operator C the coefficients C{M) are uniquely defined by C{M) = 

{n\aMC\n). 



Claim 1. Any state {ip) satisfying = 1 can be uniquely written as {ip) = exp{—C) \Q) for some 

creation operator C . 

Remark: the exponent above is defined by its Taylor series. The nilpotence of operators aj\^ 
implies that C'^ = for any k greater than the number of qubits n = \C\, so the Taylor series can be 
truncated at k = n. 

Proof: Clearly, the states {a\,^ \^)}, M (1 C, constitute the orthonormal basis of the n-qubit Hilbert 
space. Let {ip) = J2mcc'^(^)^m\^)- Equation lip) = exp (— C) \Q) is equivalent to a system of 
equations 

7/>(0) = l, C{M) = -^{M) + J2^^ Yl C{M,)---C{Mk), MCC, M^%. (7) 

k=2 ' M=A/iU...UAffc 

Here the second summation is over all partitions of M into k disjoint non-empty sets Mi, . . . , M^. 
Suppose we have already found all coefficients C{M) with \M\ < p. Then Eq. ([7]) assigns a unique 
value to all coefficients C{M) with \M\ = p+1. Thus the system Eq. ([7]) has a unique solution. □ 



2.2 Ansatz for the ground state 

Our goal is to find an eigenvector \tp) satisfying H{e) \ip) = E{e) \ip), where E{e) is the smallest 
eigenvalue of H{e). We shall use the following ansatz for \ip) (we don't care about the normalization): 

\^)=exp{-C)\n), C= J2 C{M)al. (8) 

Claim [T] asserts that the ground state can be represented in this form unless it is orthogonal to 
Since we don't require \ip) to be a normalized state, the ansatz Eq. ([H]) is meaningful only if C is a 
bounded operator. We shall define a norm of a creation operator as 

\\C\\, = maxJ2\C{M)\. (9) 

Thus the ansatz Eq. must be supplemented by a requirement that C is a creation operator with a 
finite norm ||C||i. Of course, it may happen that H has several eigenvectors of the form Eq. ([H]). One 
has to invoke some extra arguments to select an eigenvector corresponding to the smallest eigenvalue, 
see subsection 12.41 

The "physical meaning" of the norm ||C||i can be illustrated by considering a product state: 
IV') = <S>u&c 1^")' where \ipu) = |0) + a„ Obviously, \tp) = exp (-C) \n) with C = -Y^u^c^^^l- 
Accordingly, ||C||i = max^,g£ Thus one can think about ||C||i as a density of quantum fiuctua- 
tions. 

Using the identity exp (C) exp (— C) = / valid for arbitrary operator C, see [20], one can rewrite 
the Schrodinger equation H{e) \ip) = E{e) \ip) as 

exp {C)iHo)\n) + e exp iC)iV)\n) = E{e) (10) 
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Here we introduced a superoperatoi^l C such that 



C{X) = cx- xc. 

The exponent exp [C) is defined by the Taylor series. The advantage of the ansatz Eq. ([8]) is that 
we can truncate expansion of the exponent exp (C) after a few lowest orders since all higher order 
terms turn out to be identically zero. It follows from the two lemmas stated below. 

Lemma 1. Let Ci, C2 he creation operators. Then 

CiC2{Ho) = 0. (11) 

Proof: To simplify notations we shall consider operators a instead of Let u,Mi,M2 C C and 
X = aMiO'M2(|l)(lU)- By linearity, it is enough to prove that X = 0. Since the operators and 
commute, X = unless u G Mi fl M2. Then [aM2, = C'Ph = [(^Mi, '^Af2] = 0. 

□ 

Lemma 2. Let Ci, C2, ...,6*5 be creation operators. Then 

C,C2---C,iV)=0. (12) 

Proof: To simplify notations we shall consider operators a instead of a^. Let Mi, M2, . . . , M5 C C, 
{u,v) G £ and X = cim^ " " ■ (iMr,(yu,v)- By linearity it is enough to prove that X = 0. Since the 
operators qmi, ■ ■ ■ jO-Ms commute with each other, X = unless each of the subsets Mj contains at 
least one of the vertices u, v. Therefore, expanding the commutators one can represent X as a linear 
combination of 2^ terms, where each term contains at least five operators au, on the pair of qubits 
u, V. Some of these operators a are on the right of Vu^v and some of them are on the left. Thus at 
least three operators a are on the same side of Vu^v Then at least two operators a act on the same 
side of Vu^v and on the same qubit. Thus each of the 2^ terms in X contains either or aj. Thus 
X = 0. 

□ 

Combining Lemmas [TO we get the following truncations: 

exp{C)iHo) = Ho + C{Ho), (13) 

A;=0 

Here a convention C^{V) = V is adopted. 

Let us point out an analogy between the truncation effect observed above and the Lieb- Robinson 
bound [T71 [18]. The latter asserts that for any local observable acting only on a qubit u and 
for any Hamiltonian H with short-range interactions of bounded norm the time evolved observable 
Ou{t) = exp {iHt){Ou) can be approximated very well by an operator acting only on spins within 

^In our context a superoperator is a linear operator acting on the space of linear operators on Ti.. Throughout the 
paper we shall use a notation A for a superoperator A{X) = AX — XA associated with a linear operator A. 
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distance v\t\ from u, where f is a group velocity. If one takes a creation operator C for which 
the coefficients C{M) are non-zero only for subsets M of size 0(1) (an analogue of short-range 
interactions), then the "time-evolved" observable exp (C)(0„) acts only on the spins within distance 
0(1) from u (apply the same arguments as in the proof of Lemma [2]). As opposed to the Lieb- 
Robinson bound scenario, the size of a region acted on by the evolved operator does not depend on 
the norm of C (which is analogous to the evolution time ) and no approximations are involved. 

2.3 Kirkwood-Thomas equations 

Substituting Eqs. f[T51) into the Schrodinger equation Eq. flTUl) and taking into account that Hq\Q) = 
one gets 

- J2 C{M)Hoal\n) + eexp{C){V)\n) = E{e)\n). (15) 

Let us introduce eigenvalues of the unpertubed Hamiltonian Eo{M) such that 

Hoal\Q)=Eo{M)al\Q), Eo{M) = A„. (16) 

ueM 

Multiplying Eq. (fTSi) on the left by {fl\aM, M 7^ 0, and employing Eq. (IT^ one arrives at 

^^(M) = -^^5^^(n|aM0^\/)|fi), 0^MC£. (17) 

Following P], we shall refer to Eq. (fT7|) as Kirkwood-Thomas equations. Similarly, multiplying 
Eq. (fT5|) by {Q\ on the left one gets 



k=0 



It is clear that the Kirkwood-Thomas equations Eq. 0171) may have several solutions C since the 
equations do not explicitly include the eigenvalue E{e). In the worst case when neither eigenvector 
of H{e) is orthogonal to \Q) the Kirkwood-Thomas equations would have 2" solutions since any 
eigenvector could be represented in the form Eq. ([8]). We shall explain how to select the solution 
corresponding to the smallest eigenvalue in the next subsection. 

The following lemma asserts that the norm || ■ ||i has a property analogous to submultiplicativity. 
It is the main technical tool that allows one to manipulate easily with equations like Eq. ( |T71) . 

Lemma 3. Let k be any integer and Ci, . . . ,Ck be creation operators. Define a creation operator C 
such that 

C= J2 C'(M)ai, where C{M) = -^^{n\aMC^ ■ ■ ■ Ck{V)\n) . 



Then 



Eo(M) 



\ch < -^UW^Ai. (19) 
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Besides, 

k 

mC,---Ck{Vu,.)\Q)\<2'jl[\\C,\\, for any {u,v)eS. (20) 

i=i 

The proof of the lemma is presented in Appendix A. 
2.4 A lower bound on the spectral gap 

Suppose we can find some eigenvalue E'{e) of the Hamiltonian H{e) such that £"(0) = 0, E'{e) is a 
continuous function of e, and E'{e) has multiplicity 1 for |e| < ec- Then it follows immediately that 
E'{e) is the smallest eigenvalue of H{e) for all |e| < ec- Of course, the main difficulty in using this 
argument is proving non-degeneracy of an eigenvalue. The following lemma asserts that a solution 
of the Kirkwood-Thomas equations Eq. f|T71) with a sufficiently small norm ||C||i corresponds to a 
non-degenerate eigenvalue separated from the rest of the spectrum by a constant gap. 

Lemma 4. Suppose H{e) = E{e) {ip) , where {ip) = exp (— C) \Q) and C is a creation operator 
with a finite norm \\C\\i satisfying the inequality 

^ 2''dj\e\ ' iwchr .... 

k=0 

Then E{e) has multiplicity 1 and any other eigenvalue of H{e) is separated from E{e) by a gap at 
least A/2. 

Proof: Let us abbreviate H = H{e). Assume that H |0) = {E{e) + 6) |0) where \6\ < A/2 and 
the states |0) are linearly independent (the latter condition is fulfilled automatically if S ^ 0). 
We can always write |0) as 

exp (C) 10)= (22) 

MCC 

for some complex numbers B{M). Note that B{M) ^ for some non-empty set M since otherwise 
10) is proportional to {ip)- Thus we can define a creation operator B = ^0-^^fc£ -^("^) '^1/ with a 
non-zero norm > 0. Using commutativity [C, -B] = we can represent |0) as 

\^) = B\i;) + Bm\ij). 

Then the eigenvalue equations H |0) = {E{e) + 6) |0) and H \tp) = E{e) |-0) imply 

[B,H] 1^) = [B,H-Eie)I] |V^) = B{H-E{t)I) \^)-{H-E{t)I) |0) = -5 |0) = -55|^)-5i?(0) 

(23) 

Commutativity [C,B] = yields exp(C')(i?) = B. Hence, multiplying Eq. (123!) by exp (C) on the 
left one arrives at 

[B,expiC){H)] \n) + 6B\n)+6B{^) |fi) = 0. (24) 

From Lemma [1] we know that [B , exp (C) (Hq)] = [B,Ho]. Choosing any M 7^ and multiplying 
Eq. by (fi|aA/ on the left one gets 
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Here we have taken into account that BC''{V) = for k > 4, see Lemma [2l Note that condition 
\6\ < A/2 imphes a bound \Eo{M) / {Eo{M) — 6)\ < 2. Applying Lemma [3] to the operator B and 
using the triangle inequality for the norm one gets 



Bh<^-^\\Bhj2l^m\ir- 



_ k\ 

k=0 

Since ||-B||i > we can divide both sides by ||i?||i getting an inequality opposite to the one stated in 
the lemma. Thus the assumption from which we started the proof leads to a contradiction. 

□ 

Remark: Note that at e = the Hamiltonian H{e) = Hq has many degenerate eigenvalues, so one 
can certainly find two eigenvalues with separation \6\ < A/2. It might seem to be in contradiction 
with the lemma above. However at e = the condition that ||C||i is finite can not be fulfilled for 
degenerate eigenvalues, since the corresponding eigenvectors are orthogonal to 

Corollary 2. Suppose H{e) has an eigenvector \ip) = exp (— C) \Q) with an eigenvalue E\e) such 
that E\e) is a continuous function of e, E'{0) = 0, and \\C\\i < Cmax for all \e\ < e'^. Define such 
that 



k=0 



1 



Let ec = min (e^, e^'). Then for all |e| < ec 

(1) E\e) is the smallest eigenvalue of H{e) 

(2) E\e) has multiplicity 1 

(3) E'{e) is separated from the rest of the spectrum by a gap at least A/2 

Proof: (1) Indeed, Lemma H] implies that no level crossings involving the eigenvector can 
occur for |e| < ec- Since \ip) is the ground state for e = 0, it is the ground state for all |e| < ec. (2) 
and (3) follow immediately from Lemma |H 

□ 



3 Solution of the Kirkwood-Thomas equations 

In their original paper |8] Kirkwood and Thomas employed the expansion in powers of e in order 
to find the ground state. Alternative approach proposed by Datta and Kennedy in Ref. [9] and 
generahzed by Yarotsky in Ref. is to regard equation Eq. ( |T71) as a fixed point equation for a 
non-linear map on the space of creation operators. One can prove that this map is a contraction in 
the unit ball (for a properly defined metric) if e is below certain threshold value. Then one can invoke 
Brouwer fixed point theorem to argue that the unit ball contains a unique fixed point. Although the 
latter method is more elegant, we adopt the original Kirkwood-Thomas approach based on power 
series, because it naturally lends itself for getting approximation to the ground state energy with a 
controllable error. 
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3.1 Solution by formal power series 

Let us first solve Kirkwood-Thomas equation Eq. ffT7|) in terms of formal power series ignoring tlie 
convergence issue. Recall that C = J2mcc^(^) where the sum is over all non-empty sets. 
Define a series 

oo 

C{M) =Y,Cp{M)eP, 0^MC£. (25) 
p=i 

Let us agree that Co(M) = for any M. Define also 

Cp= Yl Cp{M)al, C,= J2 Cp{M)al, (26) 

SO that C = J2'^=i ^iid = J2'^=i ■ Substituting the series Eqs. ( !25]l26l) into the Kirkwood- 
Thomas equation Eq. ( fTTl) and equating the coefficients for each power of e one gets 

Ci(M) = Eo(M)-^ {Q\aMV\n), (27) 

4 1 

C,{M)=E,{M)-'Y.k\ 5Z maMC,,---C,,{V)\n), p>2. (28) 

k=l ' pi+...+Pfc=p-l 

Clearly, the equations above have a unique solution. Substituting Eqs. fl25|26l) into the formula for 
the ground state energy Eq. (fTS!) one gets 

oo 4 

E{e) = J2Epe^, E, = {n\V\n), E, = Y^- J] {n\C,,. ■ ■ C,,{Vm, p>2. (29) 

p=l k=l ' pi + ...+Pk=P~l 

Of course, formal power series do not represent an actual solution of the Kirkwood-Thomas equations 
unless we prove their convergence. 

3.2 Convergence of C-series 

We would like to prove that the series C = Yl'^=i ^p convergent with respect to the norm 

Eq. with a non-zero convergence radius. We shall need to get a lower bound on the convergence 
radius in terms of A, d, and J. Clearly, it is enough to analyze convergence of the series 

oo 

x(e) = 5Zxpe^ Xp=\\Cp\\i. (30) 
p=i 

Note that < x(|e|). 

Lemma 5. The series x(e) = J2'^=iXp^^ converges absolutely for 

H < 2.0 = — . (31) 

Besides, for any e as above one has the following bounds: 

|x(e)|<2-^^ Xp<-^ for p>l. (32) 
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Proof: Let us first get an upper bound on ||Ci||i. From Eq. (!27|) it clear that Ci(M) = unless 
M C {m, f} for some edge (m, f) G £. Let m G £ be the vertex that achieves the maximum in 
||Ci||i = max„g£ Sm9m Then the sum contains at most d + 1 sets M, namely, M = {u} 

and M = {u,v} for {u,v) G £. Therefore ||Ci||i < {d + 1)J/A < 2dJ/A, that is 

Xi < (33) 
Define a polynomial function Fp of real variables xi, . . . , according to 

P1+P2=P— 1 P1+P2+P3=P-1 pi + ...+P4=p-l 

(34) 

Applying Lemma [3] and triangle inequality to Eq. (j28l) one gets 

Xp < ^p(Xi,---,Xp-i), P>2. (35) 

To simplify notations, define constants 

2dJ 2'^dJ 

so that Xi ^ Xp ^ bFp{xi, ■ ■ ■ , Xp-i) ioi p >2. Consider the formal power series 

oo 

/i(e) = ^ /ipe^, /ii=a, iJ,p = bFp{i^i, . . . , Hp_i), p > 2. (37) 
p=i 

Since the polynomial Fp has non-negative coefficients one can prove inductively that Xp ^ A^p fo^' 
p > 1. Hence it suffices to prove that the series Eq. (157|) converges absolutely for |e| < 2eo. 

Our strategy will be to guess a function fi{e) analytic for |e| < 2eo whose Taylor series at e = 
coincides with the series Eq. fl37|) . By inspecting the recursive relation Eq. fl37|) one can easily 
convince oneself that /i(e) has to obey the following equation 

/x(e) = ae + be (^^(e) + ^fi'{e) + ^fi^e) + ^/^'(e)) ■ (38) 

We can use it to write down the inverse function 

^...N „ , . A. , 1 ..2 ,1.3 , 1 ,.4 



Simple algebra shows that 



l^?(/^)l>2 if l-^l ^ i6 = 2"''- (40) 



Thus e(/i) is analytic for < 2 . Define a set M = {/x : < 2 }. 

Claim 2. Le^ e be a complex number such that \e\ < 2eQ. Then equation e{fi) = e has a unique 
solution ji E M . 
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Proof: One can easily show that for any /ii,/i2 G M 

|Q(/^i)-Q(/W2)| < 261^1 -/i2|. (41) 

Assume e(/ii) = e(/i2) = e for some ^1,^2 € M. If e = then = fi2 = 0. Assume e 7^ 0. Then 
fj,i, /i2 7^ and 

_ /^2(Q(/il) - Q(/i2)) 

Applying the lower bound Eq. (HOj) and the upper bound Eq. (HTl) we get 

, ^ 2-i5|g(/ii)-g(/i2)| ^ 2-iH|/ii-^2| ^ i| I 

l/il - /i2| < < < -1/^1 - /i2|. 

Thus fii = fi2 and equation e(/i) = e has at most one solution /j, e M. Therefore e : M — e(M) is 
an injection. Let us prove that e(M) contains a ball of radius 2eo. Indeed, e(M) is an open set and 
G e(M). Let 7 be the boundary of M, i.e., a circle of radius 2^^^ centered at 0. Then 6(7) is the 
boundary of e(M). For any /i G 7 one has |Q(/i)| < a + 26|/i| = a + 2~^% < 2a. Thus e(M) contains 
a ball of radius 

2-15 2^^" A 
R = min |e(^)| > — — = — -— = 2eo 
f^e-y 2a aJ 

It completes the proof of the claim. 

Let K = {e : |e| < 2eo}. From Claim [2] we infer that e(/i) is an analytic bijection from the 
set e~^{K) C M to the set K. It follows from the inverse function theorem for analytic functions, 
see pro], that the inverse function fi{e) is analytic for e ^ K. Therefore the series Eq. ( l37|) converges 
absolutely for |e| < 2eo. 

The upper bound on /ip can be obtained by standard methods using Cauchy's formula: 

1 / fi{e)de 

Thus 

1 2-15 

Mol < 7 ^ niax |u(e)| < ^. 

Recall that Xp ^ A^p) so the lemma is proved. 

□ 

One can summarize the results of this subsection as follows. 



Corollary 3. Suppose \e\ < 2eQ. Then the Kirkwood- Thomas equations Eq. [F/y have a unique 
solution C defined by the power series Eq. ( fl3]) with \\C\\i < 2~^^ . 

3.3 Convergence of E'-series 

In this subsection we analyze convergence of the series E{e) = ^^e eigenvalue obtained 

from the Kirkwood-Thomas equation, see Eq. (jlj). 



14 



Lemma 6. The series E{e) = Yl^=i-^p^^ converges absolutely for 

2-i'A , , 

H < 2^0 = — . (42) 

Besides, 

\K\ < —, ^ for p > 1. (43) 

Proof: Applying Lemma [3] and triangle inequality to Eq. fl29l) one gets 

\E^\<ndJ, \Ep\<2^ndJFp{xi,---,Xp-i), P > 2, (44) 
where Xp = IIC'pHi and the polynomial Fp is defined in Eq. flMl) . Define a formal series 

oo 

e(e)=^epe^, ei = ndJ, Cp = 2'^ndJFp{xi, . . . ,Xp-i), P > 2. (45) 
By definition, \Ep\ < Cp for all p. Besides, e(e) can be expressed in terms of x(e) = Xl^i Xp^^ as 
e(e) = 2W6 (^x(e) + lx\e) + ^xHe) + ^x\e)^ + ndJ e. 

This equality can be verified by equating coefficients for each power of e. Lemma [5] implies that x(e) 
is analytic for |e| < 2eo. Therefore e(e) and E{e) are analytic for |e| < 2eo and the first statement of 
the lemma is proved. In order to get an upper bound on Cp (and thus on Ep), use Cauchy's formula: 



6p 



e{e)de 



It follows from Lemma [5] that \x{^)\ ^ 2^^^ for |e| < 2eo. Therefore 

1 1 2~-^^nA 

l^'l (2J^,SJ^'^>I - (2W(2.o)2-» + ndJ(2<„)) < 

The lemma is proved. 



□ 



Corollary 4. Suppose |e| < 2eo. T/ien t/ie series Eq. [29\) converges absolutely to the smallest 
eigenvalue of H{e). The smallest eigenvalue is non- degenerate and is separated from the rest of the 
spectrum by a gap at least A/2. 

Proof: It follows from Corollary [2l Indeed, we have already shown that the conditions of Corollary [2] 
are satisfied with = 2eo and Cmax = 2~^^, see Corollary [3l It yields e'^ < 2~^^A/{dJ). Thus 
ec = min (e^, e") = 2eo. 

□ 
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Lemma [6] allows one to estimate an error resulting from truncation of the series for the ground 
state energy at a finite order p. 

Corollary 5. Suppose \e\ < eo- Then 

V 
q=l 

Proof: Use Eq. 

□ 

Summarizing, we have proved Theorems [Tpl 

4 Linked cluster theorems 

Throughout this section we shall use a term linked cluster which refers to a subset of vertices inducing 
a connected subgraph of Q. More formally, 

Definition 2. A subset M ^ C is called a linked cluster iff for any u,v E M. there exists a sequence 
of vertices uq, ui, . . . , Ut & M such that uq = u, Ut = v and {uj, G £ for all j = 0, . . . ,t — 1. 

Definition 3. A connected size of a subset M (1 C is the minimal size of a linked cluster that 
contains all vertices of M. We shall denote a connected size of M as \M\c. 

4.1 Linked cluster expansion for the ground state 

Lemma 7. Let C{M) = Yl'^=i^pi^)^^ solution of the Kirkwood- Thomas equations obtained 

in Section\^ Then 

Cp{M) = unless \M\^ <p+l. (47) 

Proof: We shall prove the lemma by induction in p. From Eq. fl27j) one infers that Ci{M) = 
unless M C {u,v} for some edge {u,v) G S. In particular, Ci(M) = unless |M|c < 2. It proves 
the statement of the lemma for p = 1. Suppose the statement is proved for p = 1, . . . ,q — 1. From 
Eq. fl28|) one infers that Cq{M) is a linear combination of terms like 

X = Cp,{Mi) ■ ■ ■ Cp^{Mk){n\aMal^ ■ ■ ■ al^{Vu,.)\n), (48) 

where pi + . . . + pk = q — Let us figure out under what circumstances the matrix element in 
Eq. (I48p can be non-zero. 

Claim 3. Lei M, Mi, . . . , C £ be non-empty sets, iV = Mi U . . . U M^, {u, v) eS. Denote 

y = {n\aMa\^, ■ ■ ■ a\f^{Vu,v)\^) ■ 

Then y = unless the following conditions are met: 

(i) Each set Mi, . . . , M^ contains at least one of the vertices u,v. 

(ii) N\{u,v} CMC NU{u,v}. 



(46) 
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Remark: this claim is true even for M = if one adopts a convention 00 = /. 
Proof: Suppose some set Mj contains neither u nor v. Then a\j, commutes with Vu^^ as well as with 

all operators for i 7^ j. Thus y = 0. Suppose condition N\{u^v} C M is violated. Then there 
exists a set Mj and a vertex w G Mj such that w ^ u,v and w ^ M. Thus ajyj_^. contains a factor aj^ 
which commutes with all other operators involved in y. By moving a\j leftwards one can show that 
each of 2^ terms in y starts from that is, y = 0. Suppose condition M C NU{u, v} is violated, 

that is, there exists a vertex w & M such that w ^ N and w u,v. Then the operator qm contains 
a factor which commutes with all other operators involved in y. By moving rightwards one 
can show that each of 2^ terms in y tails with \ fl), that is, y = 0. 

□ 

Returning to Eq. ( HHl) we conclude that a; = unless each set Mj contains either u or/and v, and 
M C U {u,v}. Let Mj be a linked cluster of minimal size containing Mj, that is, \Mj\c = \Mj\. 
Let iV = Ml U . . . U Mfc and C = N U {u}U {v}. Then C is a linked cluster and M C C. Note that 

k k 

\C\ < \Mj\ + 2-k = 'Y\Mj\c + 2-k 
j=i i=i 

where we have taken into account that each Mj contains either u or/and v. By induction hypothesis 
we have Cp^{Mj) = unless \Mj\c < Pj + 1- Thus for any non-zero term x one has 

k 

\C\ < ^iPj + l) + 2- k = q- l + k + 2- k = q + l. 
i=i 

Thus Cq{M) = unless \M\c <q + l. □ 



4.2 Upper bound on the number of linked clusters 

The following lemma asserts that the number of linked clusters of size p containing a given vertex 
grows at most exponentially with p (if the maximal degree of the graph d is a constant). To the best 
of our knowledge, this lemma has been originally proved in Ref. [21] by Aliferis, Gottesman, and 
Preskill in the context of fault-tolerant quantum computatioijfl. 

Lemma 8. Let Np{u) be the number of linked clusters with p vertices that contain a vertex u and 
Np = maXutzc Np{u) . Then 

Np < {Ady-\ (49) 

Proof: Let Tp{u) be a set of trees with p vertices that contain a vertex u (naturally, we consider 
only those trees that are subgraphs of Q). Let Tp{u) = \Tp{u)\ be the number of such trees. For any 
tree T G Tp{u), a set of vertices of T is a linked cluster that contains u. Conversely, if M 9 m is 
a linked cluster, \M\ = p, consider a subgraph Gm induced by M. Then any spanning tree of Gm 
belongs to Tp{u). Thus Np{u) < Tp{u). 



^The authors became aware of it after completion of the present work. 
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Denote Tp = maXufzcTplu). Obviously, Ti = 1 and T2 < d. Let us prove that 

Tp<d J2 ^Pi^P- (50) 



Pl+P2=P 



where the convention Tq = is adopted. Indeed, for any edge e incident to a vertex u define a set 
Tp{u,e) that includes all trees T G Tp{u) that contain an edge e. Let Tp{u,e) = \Tp{u,e)\. Clearly, 

%){u) = LieTp{u,e), Tp{u) < y Tp{u,e) < dmaxTp{u, e). (51) 

e 

Let e = {u,v) be the edge that achieves the maximum. Note that any tree T G Tp{u,e) consists of 
the edge {u,v) and two disjoint trees Ti G Tp^{u) and T2 G Tp^{v), where Pi +P2 = Thus we have 
an upper bound 

Tp{u,e)< Tp^{u)Tp.^{v) < ^ Tp^Tp^. 

Pl+P2=P Pl+P2=P 

Substituting it into Eq. (151 p and taking the maximum over u & C we obtain Eq. (1501) . 
Define a sequence Si, S2, ■ ■ ■ such that 

Si = 1, Sp = d Sp-i^Sp2 for p > 2. (52) 

Pl+P2=P 

Clearly, Ti = Si = 1 and T2 < d = S2- It follows that Tp < Sp for all p. In order to derive 
an explicit formula for Sp define a generating function S{x) = Yl^=i '^p obeys an equation 

S{x) = dS{xY + X, which implies 

- 

Taking the derivatives one gets 



^ _l d^S 



^ p\ dxP 

It follows that 



{4dy ^ / 1 



n 



2<i(p!) J-i V 2 



a:=0 V^^V a=o 



a 



(4rf)P(p- 1)! ^ (4rf)P-i 
- 4rf(p!) - p 



Sp < ^^-^^^L/" < < (4rf)^-^ 

Summarizing, Np < Tp < Sp < {4:d)P-\ 



□ 



4.3 Linked cluster expansion for the ground state energy 

This subsection provides the necessary tools for computing spin-spin correlators. A reader interested 
only in computing the ground state energy can safely skip it. 

Let us consider a more general family of Hamiltonians for which the parameter e may be different 
on different edges. Let variable e^^v be assigned to an edge {u, v) G S. For any subset of edges ACS 
denote e[A] a collection of variables assigned to edges of A. The Hamiltonian is 

H{e[£]) = Ho+ J2 (53) 

iu,v)e£ 
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Let i?(e[£']) be the ground state energy of H{t[£]). We shall consider multivariate Taylor series for 
the function E{e[£]). 

Lemma 9. The multivariate Taylor series for E{e[£]) at the point e[£] = converges absolutely if 
\^u,v\ < Co for all {u, v) E £. 

Proof: Let Q = {e[£] : \eu,v\ < 2eo for all {u,v) G £}. Let us firstly show E{e[£]) is analytic 
function of each individual variable eu,v for e[£] G Q. Indeed, let £^ be a set of all edges except {u,v). 
Define an unperturbed Hamiltonian Hq = Hq + '^(uv)€e ^u,vVu,v and a perturbation eu,vVu,v It 
follows from Corollary H] that Hq has non-degenerate ground state and the spectral gap at least A/2. 
Applying the standard perturbation theory to a perturbed Hamiltonian Hq + eu,v Vu^v we conclude 
that E{e[£]) is analytic function of e„,^ as long as the Weyl condition He^^i, V^,t,|| < A/4 is satisfied. 
Since we assumed that \eu,v\ < 2eo, one has He^^^ < 2eoJ = 2~^'^A/d < A/4. Thus E{e[£]) is 
analytic in fl with respect to each individual variable e„^„. Repeatedly using Cauchy's formula one 
gets 

Since Hq and V are bounded operators, the absolute value can be bounded by a constant 

(maybe depending on n). The Taylor series in eu,v at the point eu,v = for any factor l/{zu,v — eu,v) 
in Eq. (15^ converges absolutely as long as |e„^„| < 2eo. Thus the Taylor series for E{e[£]) converges 
absolutely if le^^^l < eo for all {u,v) G £. 

□ 

The Taylor series for i?(e[£^]) can be uniquely written in the form 

Em) = T.\ n 'uApA{e[A]), (55) 

ACS \(u,v)eA j 

where the sum is over all subsets of edges A and p^(e[74]) is the series that involves only variables 
Eu^v pertaining to A. Clearly, the coefficients of p^(e[yl]) are functionals of interactions Vu^v with 
(m, v) E a only. The main goal of this section is to show that the expansion Eq. fl55l) involves only 
linked clusters of edges. Let us firstly define this notion. 

Definition 4. A subset of edges A C £ is called a linked cluster iff the subset of vertices induced by 
A is a linked cluster. 

Lemma 10. The series Eq. 153]) involves only linked clusters of edges A. 

Proof: Suppose A C £ is not a linked cluster of edges. Let M C £ be a set of vertices induced by 
A. Since M is not a linked cluster, it can be represented as a disjoint union M = Mi U M2, where 
Ml, M2 C £, Ml n M2 = 0, and no edge connects Mi and M2. Accordingly, A can be represented as 
a union A = Ai U A2, where Ai and A2 are the set of edges inducing Mi and M2 respectively. Let 
us choose variables e[£] such that e^^^ = unless {u,v) G A. Then it is clear that the Hamiltonian 
H{e[£]) splits into a sum of three terms acting on non-overlapping sets of qubits: 

Hie[£]) = Hi + H2 + Helse, Hj= J2 ^u\l){l\u+ ^n,vVu,., H^lse = A„|l)(l|„. 
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The ground state energy of H{e[£]) is equal to the sum of ground state energies of ifi, and Heise- 
It imphes that 

Em) = E{e[A,]) + E{e[A,]). (56) 

If we assume that pyi(e[A]) 7^ 0, when E{e[£]) would include at least one monomial including variables 
from both sets Ai, A2 which contradicts to Eq. ( l56i) . 

□ 

The following implication of Lemma [TO] will simplify computation of spin-spin correlators . 

Corollary 6. Consider a Hamiltonian H = Hq + eV, where V = J2(uv)££^u,v Let E{e) = 
Yl^=i-^p^^ series for the ground state energy of H. Suppose the interaction Vg^t depends 

on a parameter rj for some edge {s,t) G S. Then a derivative 



K. - 



ri=0 



can be computed by setting Vu^v = for all edges {u,v) having distance p+ 1 or greater from the edge 
{s,t). 

Proof: Indeed, Ep can be obtained from Eq. (l55l) by setting eu,v = e on all edges, restricting the 
sum to linked clusters A of size at most p and collecting all monomials of total degree p. Clusters A 
that do not contain the edge (s, t) will not contribute to Kp. Clusters A that contain the edge (s, t) 
cannot contain any edge (m, v) having distance p + 1 or greater from the edge (s, t). 

□ 

5 Computational algorithms 

In this section we describe an algorithm that takes as input a description of the Hamiltonians Hq, V 
and an integer p. The algorithm returns a list of coefficients Ei, Ep in the series for the ground 
state energy E(e) = Yl'^=i -^p The running time of the algorithm is nexp (0(p)). In Section [5751 we 
describe a generalization of the algorithm that allows one to compute spin-spin correlation functions. 



5.1 Computing the coefficients Cp{M) 

The first part of the algorithm is to compute the coefficients Cq(M), ^ ^ M C sequentially for 
q = 1, . . . ,p using the solution of the Kirkwood-Thomas equations, see Eqs. fl27ll28p . This gives an 
approximate description of the ground state. 

We shall store triples {M,q,Cq{M)) in n bins (memory registers) Bu labeled by vertices of the 
graph u & C. Once a coefficient Cg{M) is computed, the triple {M,q,Cg{M)) is placed into every 
bin Bu for which u G M. From Lemma [7] we learn that Cq{M) = unless M is a subset of some 
linked cluster M of size at most q + 1. According to Lemma [HI the number of linked clusters M 
of size g + 1 containing vertex u is bounded by exp (0(g)), where the coefficient in the exponent 
depends only on d. Each linked cluster of size q + 1 containing vertex u has 2"^ subsets containing 
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vertex u. Thus we can bound the number of entries in the bin Bu at the moment when all coefficients 
Ci(M), . . . , Cp{M) have been computed as < J2q=i 2''exp (0(g)) = exp (0(p)). 

Suppose we have already computed all non-zero coefficients Ci(M), . . . , Cg_i(M), M C £. The 
next step is to compute coefficients Cg{M) for all sets M C C satisfying the condition of Lemma [TJ 
that is |M|c < q + 1. Expanding Eq. (1251) one gets 

{u,v)££ k=l ' pi+...+Pfc=g-l Mi,...,MfcC£ 

(57) 

Note that the right hand side of this equation involves only coefficients Cp.{Mj) that have been 
already computed. For simplicity let us assume that computation of any term in Eq. fl571) requires 
one unit of tim^. Denote 

X = {niaMttlj^ ■ ■ ■a\^^{Vu,v)\^)- (58) 

Recall, see Claim [3l that x = unless the following conditions are met: 

(i) Each set Mi, . . . , contains at least one of the vertices u, v. 

(ii) N\{u, v} C M C N U {u, v}, where N = M-^U . . .U Mk. 

The property (i) implies that for a fixed edge {u, v) we can restrict the three rightmost sums in 
Eq. (1571) by taking triples {Mj,pj, Cp.{Mj)) either from the bin or from the bin B^. Therefore for 
a fixed {u,v), the overall number of non-zero terms in the three rightmost sums in Eq. f l57p can be 
bounded by (|5,| + \B,\)'^ < (|5„| + \B,\)^ = exp (0(g)). 

We shall now prove that only a small number of edges {u,v) can give a non-zero contribution to 
Cq{M). Indeed, there are two cases: (1) M C {u,v}] (2) There exists w E M such that w ^ {u,v}. 
Clearly only 0(1) edges {u,v) can lead to the case (1), so let us focus on the case (2). Consider 
any term x as in Eq. ( l58l) . Properties (i),(ii) above imply that x = unless there exists a set Mj 
such that w G Mj and one of the vertices u,v belongs to Mj. Without loss of generality, w,u E Mj. 
Lemma [7] implies that Cp.{Mj) = unless \Mj\c < + 1 < g. Therefore the distance between u 
and w is at most g. Taking into account that \M\ < \M\c < g + 1, we can bound the number of 
edges {u,v) that can give a non-zero contribution to Cg{M) by |M|(i''+-'^ < {q + l)d'^'^^ = exp (0(g)). 
Summarizing, the overall number of non-zero terms in Eq. (!57|) is exp (0(g)). 

In order to compute all non-zero coefficients Cq{M) we will have to repeat the procedure above 
for each subset M satisfying the condition of Lemma [TJ that is \M\c < g + 1. By Lemma [HI the 
number of such sets is n exp (0(g)). Summarizing, the overall time one needs to compute all the 
coefficients C'i(M), . . . , Cp{M) is nexp (0(p)). 

5.2 Computing the ground state energy 

The final step of the algorithm is to compute the coefficients Ei, . . . ,Ep using Eq. (l29l) . This equation 
can be expanded as 

4 1 

^P= E E Cp,iM,)---Cp,{M,){Q\aM,---aM,{V.M. (59) 

iu,v)e£ k=l ' pi + ...+Pfc=p-l Afi,...,MfeC£ 

^ This assumption might seem mijustified, because the precision up to which the coefScient Cq{M) must be 
computed depends upon S. However, taking into account these subtleties will lead to an additional overhead 
poZy(logn, log(5^^) which can be neglected since the algorithm has running time poly{n, S~^). 
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Expanding the commutators in the matrix element one gets 2^ terms. However, only the term in 
which Vu^v is the leftmost operator gives a non-zero contribution. Thus 

(Sl\aM, ■ ■ ■aM,(K,.)|^^) = (-1)' W^u^vO^M, ■ ■ -aUl^)- 

In particular, we can restrict the summation over Mi, . . . , by subsets of {m, f } only. There are 
only 3 such subsets: {u}, {f}, and {u, f}. This observation implies that the number of terms in the 
rightmost sum in Eq. fl59|) is 0(1). Since there are Oij?^ partitions pi + . . . + = P ~ 1? ^ < 4, the 
overall number of terms in Eq. fl59l) is O(np^). 

Combining the results of Subsections 15. 1|5. 21 we conclude that the overall time needed to compute 
the coefficients Ei, . . . ,Ep scales as ?7,exp (0(p)). In the above analysis we did not keep track of the 
coefficients in the exponents 0{p). If one computes the exact coefficient, it yields the overall running 
time 77,2i^+^'°s('^), where log stands for the base two logarithm. Accordingly, the running time as a 
function of n and 6 scales as 

For example, implementing the algorithm on a 2D square lattice {d = 4) would require a running 
time T{n,6) ~ n(n(5~^)^^, which is certainly not practical. Note however, that the power of n6~^ 
depends upon the ratio \e\/R, where R is the convergence radius of the series Eq. (jl]). The power 
15 + 6 log (d) corresponds to the most pessimistic scenario R = 2eo (the best lower bound on the 
convergence radius that we can prove) and |e| = eo. 



5.3 Computing spin-spin correlation functions 

Let s, t G £ be any pair of vertices. It may or may not be an edge of the graph Q. Let us add {s, t) 
to the set of edges S (by creating a double edge between s and t if necessary). The modified graph 
has maximal degree d* = d + 1. Let Os,t be a Hermitian operator acting non-trivially only on qubits 
s,t. We shall assume that ||Os,f|| < J. The quantity we are interested in is the expectation value 

where {ip) is the ground state of H{e) = Ho + eV. Our goal is to compute K with a specified precision 
6. To this end let us define a Hamiltonian 



H{e,r])=Ho + eV + r]Os,t. (60) 
Let E{e,ri) be the smallest eigenvalue of H{e,ri). As we know from Lemma [HI the Taylor series 

oo 

Eie,v)= Y,E,,,e^r^^ (61) 

p,q=0 

converges absolutely for |e|, \ri\ < eg, where 

2-18 A 2-18 A 



d*J {d+1) J' 
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The Hellman-Feynman theorem asserts that 



K 



drj 



ri=0 



p=0 



Our algorithm will get an approximation to K by computing a truncation of series in Eq. 
following lemma provides a bound on the error resulting from the truncation. 

Lemma 11. Suppose \e\ < eQ/{2d). Then 



(62) 



The 



q=0 



< 2 



-16-p 



Jd{d+l). 



Proof: Let us firstly prove that 



)-16 M+l 



IK 



< 



Indeed, use Cauchy's formula: 



(2'Kl) 



E{e, rj) de drj 



-1 



(63) 



(64) 



(65) 



From Lemma [6] we infer that |-E'(e, r^)! < 2^^^nA. However, we would like to have an upper bound 
independent of n. To this end we employ Corollary [6] according to which Ep^i can be computed 
by restricting the Hamiltonian on [d + l)-neighborhood of the edge {s,t). The number of spins 
in this neighborhood is at most n* = dP~^^. Therefore \E{e,7])\ < 2~^^ d^^^A. Substituting this 
bound into Eq. fl65l) one gets Eq. (IMl) . Finally, using the condition |e| < eQ/{2d) we bound the sum 
E^p+il^.,i|e^ as in Eq. ([63]). 

□ 

Lemma [11] shows that in order to compute K with an absolute error 6 it is enough to compute 
the coefficients Eq^i, £'i i, . . . , Ep^i in the series Eq. ( 1611) with p = log (5~^) + 0(1). 

Computation of the coefficients Ep^i requires only minor modifications of the algorithm described 
in Sections I5.1f5.2[ Indeed, consider a function E{e,ri) = E{e,eri). Using the series Eq. one gets 



p+q= 



In particular. 



dEp+i{r]) 



drj 



(66) 



»7=0 



On the other hand, i^(e, rj) is the ground state energy of a Hamiltonian Hq + e {V + rj Og^t)- Thus we 
can compute the coefficients £'1(77), ... , Ep^i{ri) using the already available algorithm for the ground 



Ep^i can be 



state energy. Moreover, from Corollary [6] we know that the coefficients i?o,i, • 
computed by restricting the Hamiltonian to the {p + l)-neighborhood of the edge (s,t). Thus we 
can apply Theorem [3] with n replaced by n* = dP~^^, obtaining an algorithm with a running time 
exp(0(p)) for computing Ei{ri) , . . . , Ep^i{ri) . In fact, at every step of this algorithm we have to 
retain only the terms independent of rj and the terms linear in t], see Eq. ( |66|) . Since we have chosen 
p = log + 0(1), the running time of the algorithm is poly{5^^). 
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6 Discussion and open problems 



We have proved that the ground state properties of a spin Hamiltonian with sufficiently weak inter- 
actions between qubits can be computed efficiently. We hope that this result could be generalized 
in several different directions. Firstly, one could try to consider more general class of unperturbed 
Hamiltonians Hq, for example, classical Ising-like Hamiltonians. In addition, one could consider sys- 
tems of fermionic modes rather than spins. Secondly, one could investigate possible generalizations 
of the Kirkwood-Thomas ansatz to the case of degenerate ground state. In this case the ansatz 
should be constructed for an effective Hamiltonian acting on a low-energy subspace rather than for 
the ground state. Results of this kind could provide a rigorous basis for perturbative derivations 
of low-energy effective Hamiltonians, for example the mapping from the half-filled Hubbard model 
to the Heisenberg model. Thirdly, one could try to get a stronger lower bound on the convergence 
radius R of the series E{e) = Xl^i-^p^^- note that a stronger lower bound R > A/dJ can be 
easily obtained for classical Hamiltonians, when all interactions Vu^v are diagonal in the |0), |1) basis. 
Therefore, one could speculate that in the quantum case R should be close to A/dJ. 
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8 Appendix A 



In this section we prove Lemma [3l By definition of the norm, ||C||i = max^g^Fu, where 

Y^=J2Eo{M)-'\{n\aMC,---Ck{V)\Q) . 

Applying the triangle inequality one can bound F„ as 

MBu {v,w)gS Mi,...,Mk 

Here the last sum is over all non-empty subsets Mi, . . . , Mk C C. Claim [3] allows one to restrict the 
summation in Eq. fl67|) only by tuples (M, Mi, . . . ,Mk,v,w) satisfying conditions (i),(ii). We shall 
partition into A; + 1 (possibly overlapping) sums that will be dealt with separately. We define 
Xu \ j = 1, . . . , A; as a sum of all terms in Eq. (1671) for which u G Mj. We define Xu^ as a sum of all 
terms in Eq. (1671) for which u G {v,w}. In other words, 

X(^-)= ^Eo(M)-i J2 J2 XM,H|(^l«M<---aL^(V;,J|l^)| |C'i(Mi)|---|C',(M,)|, 

A/9U {v,w)e£ Mi,...,Mk 



24 



where XMj is the characteristic functioE0 of Mj and 

MBu v:{u,v)eeMi,...,Mk 

Condition (ii) in Claim [3] imphes that u G M C U {v,w}, so that each term in Eq. fl67j) appears 
at least one time in the sums Xu \ • • • , Xu^\ hence 

k 

Xu<Y^u^- (68) 

j=0 

Upper bound on X^^\ 1 < j < k: The property (i) in Claim [3] implies that at least one end-point 
of the edge {v,w) belongs to Mj. W.l.o.g. v G Mj. Then property (ii) implies Mj C M U {w}, 
so that \Mj\ < 2\M\ (recall that M is a non-empty set because u G M). It gives us a bound 
Eo{M) > A\M\ > {A/2)\Mj\. Note also that for any fixed Mi,...,Mk and v,w there exist at 
most four sets M satisfying condition (ii) of Claim [3] (take N and add/subtract vertices v and w). 
Therefore 

{v,w)e€ Mi,...,Mk ^ 

Now we can bound the matrix element by 2^ J and add a restriction Mjn{f , w} 7^ to the summations 
over sets Mj, i ^ j, see Claim |3l property (i). Taking into account that 

J2 \C^{M,)\ < J2 \C^{M,)\ + J2 \C^{M,)\ < 2\\Q\\i (69) 

we arrive to 

o2fc+2 T 1 

xl^^<^^X{\mi E Y.^M,{u)xMMr^\\CAM,)\. 

Changing the order of summations and bounding the sum over {v,w) by (i|Mj| one gets 

n2k+2j J 92fc+2 J j 

Finally, Lemma [2] implies that it is enough to consider A; < 4, so that 

k 912 J 7- 

E-^i-'S^IIllQIk. (70) 

Upper bound on X^^'^: Claim [3] implies that for any fixed (Mi, . . . , Mfc, v) there exist at most four 
sets M satisfying {ii). Using a bound Eq{M) > A we arrive to 

i);(n,i))e£Mi,...,A4 



'That is (u) = 1 if u e Afj and (") = otherwise. 
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Claim [3l property (i) allows us to bound the matrix element by 2^ J and add a restriction Miflj-u, f } 7^ 
to the summations over sets Mj. Using Eq. (1691) we arrive to 

92fc+2 J ^ nlOi T 

j = l t):{M,i))e£' j = l 

Combining Eqs. fl68|70ll71l) we prove the upper bound Eq. f|T9l) . 

The second bound Eq. fl20|) of Lemma [3] is much easier to prove. Applying the triangle inequality 
one gets 

■ ■■Ck{yu,.)\n)\ < K^I«U ■ ■ ■44(K,.)|i^)||Ci(Mi)| • • ■ \Ck{Mk)\. 

Mi,...,Mk 

Clearly the matrix element is zero unless Mj C {u,v} for all j. Expanding the commutators one 
gets 2^ terms, but only the term in which all creation operators stand on the right of Vu^v gives a 
non-zero contribution. Taking into account that 

^ |Q(M,)|<2||C,||i, 

Mj : MjC-{u,v} 

one arrives at 

k k 

m\C, ■ ■■CkiV^^.mi < 2^=7 J] < 2Vn 

i=i j=i 
where we have applied Lemma [2] to argue that k < 4. 
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